package br.edu.ufcg.msnlab2.grupo05.splineCubic;

import java.util.ArrayList;

import br.edu.ufcg.msnlab2.misc.Point;
import br.edu.ufcg.msnlab2.misc.PointImpl;


public class SplineCubicFuncao {

	private int nPoints = 0; // no. of tabulated points
	private int nPointsOriginal = 0; // no. of tabulated points after any
										// deletions of identical points
	private double[] y = null; // y=f(x) tabulated function
	private double[] x = null; // x in tabulated function f(x)
	private int[] newAndOldIndices; // record of indices on ordering x into
									// ascending order
	private double xMin = Double.NaN; // minimum x value
	private double xMax = Double.NaN; // maximum x value
	private double range = Double.NaN; // xMax - xMin
	private double[] d2ydx2 = null; // second derivatives of y
	private double yp1 = Double.NaN; // first derivative at point one
	// default value = NaN (natural spline)
	private double ypn = Double.NaN; // first derivative at point n
	// default value = NaN (natural spline)
	private boolean derivCalculated = false; // = true when the derivatives have
												// been calculated
	private String subMatrixIndices = " "; // String of indices of the
											// submatrices that have called
											// CubicSpline from higher order
											// interpolation

	private boolean checkPoints = false; // = true when points checked for
											// identical values
	private boolean averageIdenticalAbscissae = false; // if true: the the
														// ordinate values for
														// identical abscissae
														// are averaged
	// if false: the abscissae values are separated by 0.001 of the total
	// abscissae range;
	private static double potentialRoundingError = 5e-15; // potential rounding
															// error used in
															// checking wheter a
															// value lies within
															// the interpolation
															// bounds
	private static boolean roundingCheck = true; // = true: points outside the
													// interpolation bounds by
													// less than the potential
													// rounding error rounded to
													// the bounds limit

	// Constructors
	// Constructor with data arrays initialised to arrays x and y
	public SplineCubicFuncao(double[] x, double[] y) {
		this.nPoints = x.length;
		this.nPointsOriginal = this.nPoints;
		if (this.nPoints != y.length)
			throw new IllegalArgumentException(
					"Arrays x and y are of different length " + this.nPoints
							+ " " + y.length);
		if (this.nPoints < 3)
			throw new IllegalArgumentException(
					"A minimum of three data points is needed");
		this.x = new double[nPoints];
		this.y = new double[nPoints];
		this.d2ydx2 = new double[nPoints];
		for (int i = 0; i < this.nPoints; i++) {
			this.x[i] = x[i];
			this.y[i] = y[i];
		}
		orderPoints();
	}



	// Sort points into an ascending abscissa order
	public void orderPoints() {
		double[] dummy = new double[nPoints];
		this.newAndOldIndices = new int[nPoints];
		// Sort x into ascending order storing indices changes
		Fmath.selectionSort(this.x, dummy, this.newAndOldIndices);
		// Sort x into ascending order and make y match the new order storing
		// both new x and new y
		Fmath.selectionSort(this.x, this.y, this.x, this.y);

		// Minimum and maximum values and range
		this.xMin = Fmath.minimum(this.x);
		this.xMax = Fmath.maximum(this.x);
		range = xMax - xMin;
	}

	// Checks for and removes all but one of identical points
	// Checks and appropriately handles identical abscissae with differing
	// ordinates
	public void checkForIdenticalPoints() {
		int nP = this.nPoints;
		boolean test1 = true;
		int ii = 0;
		while (test1) {
			boolean test2 = true;
			int jj = ii + 1;
			while (test2) {
				if (this.x[ii] == this.x[jj]) {
					if (this.y[ii] == this.y[jj]) {
						System.out.print(subMatrixIndices
								+ "CubicSpline: Two identical points, "
								+ this.x[ii] + ", " + this.y[ii]);
						System.out.println(", in data array at indices "
								+ this.newAndOldIndices[ii] + " and "
								+ this.newAndOldIndices[jj]
								+ ", latter point removed");

						for (int i = jj; i < nP - 1; i++) {
							this.x[i] = this.x[i + 1];
							this.y[i] = this.y[i + 1];
							this.newAndOldIndices[i - 1] = this.newAndOldIndices[i];
						}
						nP--;
						for (int i = nP; i < this.nPoints; i++) {
							this.x[i] = Double.NaN;
							this.y[i] = Double.NaN;
							this.newAndOldIndices[i - 1] = -1000;
						}
					} else {
						if (this.averageIdenticalAbscissae == true) {
							System.out
									.print(subMatrixIndices
											+ "CubicSpline: Two identical points on the absicca (x-axis) with different ordinate (y-axis) values, "
											+ x[ii] + ": " + y[ii] + ", "
											+ y[jj]);
							System.out
									.println(", average of the ordinates taken");
							this.y[ii] = (this.y[ii] + this.y[jj]) / 2.0D;
							for (int i = jj; i < nP - 1; i++) {
								this.x[i] = this.x[i + 1];
								this.y[i] = this.y[i + 1];
								this.newAndOldIndices[i - 1] = this.newAndOldIndices[i];
							}
							nP--;
							for (int i = nP; i < this.nPoints; i++) {
								this.x[i] = Double.NaN;
								this.y[i] = Double.NaN;
								this.newAndOldIndices[i - 1] = -1000;
							}
						} else {
							double sepn = range * 0.0005D;
							System.out
									.print(subMatrixIndices
											+ "CubicSpline: Two identical points on the absicca (x-axis) with different ordinate (y-axis) values, "
											+ x[ii] + ": " + y[ii] + ", "
											+ y[jj]);
							boolean check = false;
							if (ii == 0) {
								if (x[2] - x[1] <= sepn)
									sepn = (x[2] - x[1]) / 2.0D;
								if (this.y[0] > this.y[1]) {
									if (this.y[1] > this.y[2]) {
										check = stay(ii, jj, sepn);
									} else {
										check = swap(ii, jj, sepn);
									}
								} else {
									if (this.y[2] <= this.y[1]) {
										check = swap(ii, jj, sepn);
									} else {
										check = stay(ii, jj, sepn);
									}
								}
							}
							if (jj == nP - 1) {
								if (x[nP - 2] - x[nP - 3] <= sepn)
									sepn = (x[nP - 2] - x[nP - 3]) / 2.0D;
								if (this.y[ii] <= this.y[jj]) {
									if (this.y[ii - 1] <= this.y[ii]) {
										check = stay(ii, jj, sepn);
									} else {
										check = swap(ii, jj, sepn);
									}
								} else {
									if (this.y[ii - 1] <= this.y[ii]) {
										check = swap(ii, jj, sepn);
									} else {
										check = stay(ii, jj, sepn);
									}
								}
							}
							if (ii != 0 && jj != nP - 1) {
								if (x[ii] - x[ii - 1] <= sepn)
									sepn = (x[ii] - x[ii - 1]) / 2;
								if (x[jj + 1] - x[jj] <= sepn)
									sepn = (x[jj + 1] - x[jj]) / 2;
								if (this.y[ii] > this.y[ii - 1]) {
									if (this.y[jj] > this.y[ii]) {
										if (this.y[jj] > this.y[jj + 1]) {
											if (this.y[ii - 1] <= this.y[jj + 1]) {
												check = stay(ii, jj, sepn);
											} else {
												check = swap(ii, jj, sepn);
											}
										} else {
											check = stay(ii, jj, sepn);
										}
									} else {
										if (this.y[jj + 1] > this.y[jj]) {
											if (this.y[jj + 1] > this.y[ii - 1]
													&& this.y[jj + 1] > this.y[ii - 1]) {
												check = stay(ii, jj, sepn);
											}
										} else {
											check = swap(ii, jj, sepn);
										}
									}
								} else {
									if (this.y[jj] > this.y[ii]) {
										if (this.y[jj + 1] > this.y[jj]) {
											check = stay(ii, jj, sepn);
										}
									} else {
										if (this.y[jj + 1] > this.y[ii - 1]) {
											check = stay(ii, jj, sepn);
										} else {
											check = swap(ii, jj, sepn);
										}
									}
								}
							}

							if (check == false) {
								check = stay(ii, jj, sepn);
							}
							System.out
									.println(", the two abscissae have been separated by a distance "
											+ sepn);
							jj++;
						}
					}
					if ((nP - 1) == ii)
						test2 = false;
				} else {
					jj++;
				}
				if (jj >= nP)
					test2 = false;
			}
			ii++;
			if (ii >= nP - 1)
				test1 = false;
		}
		this.nPoints = nP;
		if (this.nPoints < 3)
			throw new IllegalArgumentException(
					"Removal of duplicate points has reduced the number of points to less than the required minimum of three data points");

		this.checkPoints = true;
	}

	// Swap method for checkForIdenticalPoints procedure
	private boolean swap(int ii, int jj, double sepn) {
		this.x[ii] += sepn;
		this.x[jj] -= sepn;
		double hold = this.x[ii];
		this.x[ii] = this.x[jj];
		this.x[jj] = hold;
		hold = this.y[ii];
		this.y[ii] = this.y[jj];
		this.y[jj] = hold;
		return true;
	}

	// Stay method for checkForIdenticalPoints procedure
	private boolean stay(int ii, int jj, double sepn) {
		this.x[ii] -= sepn;
		this.x[jj] += sepn;
		return true;
	}



	// Calculates the second derivatives of the tabulated function
	// for use by the cubic spline interpolation method (.interpolate)
	// This method follows the procedure in Numerical Methods C language
	// procedure for calculating second derivatives
	public void calcDeriv() {
		double p = 0.0D, qn = 0.0D, sig = 0.0D, un = 0.0D;
		double[] u = new double[nPoints];

		if (Double.isNaN(this.yp1)) {
			d2ydx2[0] = u[0] = 0.0;
		} else {
			this.d2ydx2[0] = -0.5;
			u[0] = (3.0 / (this.x[1] - this.x[0]))
					* ((this.y[1] - this.y[0]) / (this.x[1] - this.x[0]) - this.yp1);
		}

		for (int i = 1; i <= this.nPoints - 2; i++) {
			sig = (this.x[i] - this.x[i - 1]) / (this.x[i + 1] - this.x[i - 1]);
			p = sig * this.d2ydx2[i - 1] + 2.0;
			this.d2ydx2[i] = (sig - 1.0) / p;
			u[i] = (this.y[i + 1] - this.y[i]) / (this.x[i + 1] - this.x[i])
					- (this.y[i] - this.y[i - 1]) / (this.x[i] - this.x[i - 1]);
			u[i] = (6.0 * u[i] / (this.x[i + 1] - this.x[i - 1]) - sig
					* u[i - 1])
					/ p;
		}

		if (Double.isNaN(this.ypn)) {
			qn = un = 0.0;
		} else {
			qn = 0.5;
			un = (3.0 / (this.x[nPoints - 1] - this.x[this.nPoints - 2]))
					* (this.ypn - (this.y[this.nPoints - 1] - this.y[this.nPoints - 2])
							/ (this.x[this.nPoints - 1] - x[this.nPoints - 2]));
		}

		this.d2ydx2[this.nPoints - 1] = (un - qn * u[this.nPoints - 2])
				/ (qn * this.d2ydx2[this.nPoints - 2] + 1.0);
		for (int k = this.nPoints - 2; k >= 0; k--) {
			this.d2ydx2[k] = this.d2ydx2[k] * this.d2ydx2[k + 1] + u[k];
		}
		this.derivCalculated = true;
	}

	// INTERPOLATE
	// Returns an interpolated value of y for a value of x from a tabulated
	// function y=f(x)
	// after the data has been entered via a constructor.
	// The derivatives are calculated, bt calcDeriv(), on the first call to this
	// method ands are
	// then stored for use on all subsequent calls
	public double interpolate(double xx) {

		if (!this.checkPoints)
			this.checkForIdenticalPoints();
		// Check for violation of interpolation bounds
		if (xx < this.x[0]) {
			// if violation is less than potntial rounding error - amend to lie
			// with bounds
			if (SplineCubicFuncao.roundingCheck
					&& Math.abs(x[0] - xx) <= Math.pow(10, Math.floor(Math
							.log10(Math.abs(this.x[0]))))
							* SplineCubicFuncao.potentialRoundingError) {
				xx = x[0];
			} else {
				throw new IllegalArgumentException("x (" + xx
						+ ") is outside the range of data points (" + x[0]
						+ " to " + x[this.nPoints - 1] + ")");
			}
		}
		if (xx > this.x[this.nPoints - 1]) {
			if (SplineCubicFuncao.roundingCheck
					&& Math.abs(xx - this.x[this.nPoints - 1]) <= Math.pow(10,
							Math.floor(Math.log10(Math
									.abs(this.x[this.nPoints - 1]))))
							* SplineCubicFuncao.potentialRoundingError) {
				xx = this.x[this.nPoints - 1];
			} else {
				throw new IllegalArgumentException("x (" + xx
						+ ") is outside the range of data points (" + x[0]
						+ " to " + x[this.nPoints - 1] + ")");
			}
		}

		if (!this.derivCalculated)
			this.calcDeriv();

		double h = 0.0D, b = 0.0D, a = 0.0D, yy = 0.0D;
		int k = 0;
		int klo = 0;
		int khi = this.nPoints - 1;
		while (khi - klo > 1) {
			k = (khi + klo) >> 1;
			if (this.x[k] > xx) {
				khi = k;
			} else {
				klo = k;
			}
		}
		h = this.x[khi] - this.x[klo];

		if (h == 0.0) {
			throw new IllegalArgumentException(
					"Two values of x are identical: point " + klo + " ("
							+ this.x[klo] + ") and point " + khi + " ("
							+ this.x[khi] + ")");
		} else {
			a = (this.x[khi] - xx) / h;
			b = (xx - this.x[klo]) / h;
			yy = a
					* this.y[klo]
					+ b
					* this.y[khi]
					+ ((a * a * a - a) * this.d2ydx2[klo] + (b * b * b - b)
							* this.d2ydx2[khi]) * (h * h) / 6.0;
		}
		return yy;
	}

	
	public String returnAllCubicPoly(){
        String answer ="";
        for (int i = 0; i < x.length-1; i++){
            answer += getCubicPoly(x[i], x[i+1], y[i], y[i+1])+"\n";
        }
        return answer;
    }
	
	
	public ArrayList<String> returnListAllCubicPoly(){
        String answer ="";
        ArrayList<String> funcoes = new ArrayList<String>();
        for (int i = 0; i < x.length-1; i++){
            funcoes.add(getCubicPoly(x[i], x[i+1], y[i], y[i+1]));
        }
        return funcoes;
    }

    
    private String getCubicPoly(double x0, double x1, double y0, double y1) {

        double alfa = 0;
        double beta = 0;
        double gama = 0;
        double h = x1- x0;
        alfa = ((y1-y0)/ h)- ((y1* h)/6 ) - ((y0* h)/ 3) ;
        beta = y0 /2;
        gama = (y1 - y0)/ (6* h);
       
        return "For x between "+ x0+" and "+x1+
        "-> pi(x)= " + y0+ " + ("+ alfa +")(x - (" +x0+ ")) + (" +beta+")(x - ("+ x0+"))^2 + ("+ gama +")(x - ("+x0+"))^3";

    }
    
    public ArrayList<Point> getPontos() {
    	ArrayList<Point> points = new ArrayList<Point>();
		for (int i = 0; i < x.length; i++) {
			points.add(new PointImpl(this.x[i], this.y[i]));
		}
		return points;
	}
}
